rm(list = ls())

library(ggthemes)
library(lubridate)
library(estimatr)
library(texreg)
library(tidyverse)

#Data should be obtained from harvard dataverse: 
#https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/0ONBF6
df = readstata13::read.dta13("data/Kao Revkin Retribution or Reconciliation AJPS 2021 Dataset.dta") 


# DATA PREP ----------------------------------------------------------------
df = df %>% 
  mutate(Act = recode(Act, `1` = "IS Fighter",
                      `2` = "Cook for Fighters",
                      `3` = "Married a Fighter",
                      `4` = "Janitor at Municipality",
                      `5` = "Paid Taxes"),
         Act = factor(Act, levels = c("IS Fighter", "Cook for Fighters",
                                      "Married a Fighter", "Janitor at Municipality",
                                      "Paid Taxes"
                                      )),
         Woman = recode(Woman, `0` = "Male", `1` = "Female"),
         Woman = factor(Woman, levels = c("Male", "Female")),
         Youth0 = recode(Youth, `0` = "Elder", `1` = "Youth"),
         Youth0 = factor(Youth0, levels = c("Elder", "Youth")),
         Tribal_Member = recode(Tribal_Member, `0` = "Not Tribal Member", `1` = "Tribal Member"),
         Tribal_Member = factor(Tribal_Member, levels = c("Not Tribal Member", "Tribal Member")))




# Figure 10A----------------------------------------------------------------

df = df %>% mutate(Volition = recode(Volition, `0` = "Involuntary",
                                     `1` = "Voluntary"),
                   Volition = factor(Volition, levels = c("Involuntary", 
                                                          "Voluntary")),
                   Act = fct_rev(Act))

mod2 = lm_robust(Punish ~ (Act + Woman + Youth + Tribal_Member)*Volition, clusters = R_ID, data = df)

library(sjPlot)
fig10a = plot_model(mod2, type = "pred", terms = c("Act", "Volition")) +
  theme_bw() + 
  coord_flip(ylim = c(1, 5)) + 
  theme(legend.position = "bottom") + 
  scale_color_manual(values = c("black", "gray54")) + 
  labs(y = "Likelihood of Punishment (Scale: 1-5)", x = "", title = "") + 
  theme(legend.title = element_blank(), axis.title = element_text(size = 9))
fig10a



# Figure 10b----------------------------------------------------------------
mod3 = lm_robust(Forgive_All ~ (Act + Woman + Youth + Tribal_Member)*Volition, clusters = R_ID, data = df)

fig10b = plot_model(mod3, type = "pred", terms = c("Act", "Volition")) +
  theme_bw() + 
  coord_flip(ylim = c(0, .9)) + 
  theme(legend.position = "bottom") + 
  scale_color_manual(values = c("black", "gray54")) + 
  labs(y = "Likelihood of Forgiveness (Scale: 0-1)", x = "", title = "") + 
  theme(legend.title = element_blank(), axis.title = element_text(size = 9)) + 
  scale_y_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8))
fig10b



library(patchwork)
fig10_complete = fig10a + fig10b + 
  plot_annotation(title = "Perceptions of Voluntariness Interacted with Enemy Act") + 
  theme(plot.title = element_text(hjust = 0.5))
fig10_complete

ggsave("figs/irq_fig10.pdf", fig10_complete, height = 5, width = 10)

# Figure 9 ---------------------------------------------------------------
df$Victimization0 = df$Victimization
df = df %>% 
  mutate(Act = factor(Act, levels = c("IS Fighter", "Cook for Fighters", "Married a Fighter",
                                      "Janitor at Municipality", "Paid Taxes")))
regb5_1 = lm_robust(Punish ~ Victimization0 * (Act + Woman + Youth0 + Tribal_Member), clusters = R_ID, data = df)
regb5_2 = lm_robust(Forgive_All ~ Victimization0 * (Act + Woman + Youth0 + Tribal_Member), clusters = R_ID, data = df)

fig9a = regb5_1 %>% 
  tidy %>% 
  filter(!str_detect(term, "Intercept")) %>% 
  mutate(term = term %>% str_remove("Act") %>% 
           str_remove("Woman") %>% str_remove("Youth0") %>% 
           str_remove("Tribal_Member") %>% str_remove("Victimization0"),
         term = recode(term, "Victimization:Cook for Fighters" = "Cook X Victimization",
                       "Victimization:Married a Fighter" = "Married X Victimization",
                       "Victimization:Janitor at Municipality" = "Janitor X Victimization",
                       "Victimization:Paid Taxes" = "Paid Taxes X Victimization",
                       "Victimization:Female" = "Female X Victimization",
                       "Victimization:Youth" = "Youth X Victimization",
                       "Victimization:Tribal Member" = "Tribal Member X Victimization"
                       )) %>% 
  mutate(term = factor(term, levels = c("Cook for Fighters", "Married a Fighter",
                                        "Janitor at Municipality", "Paid Taxes", "Victimization",
                                        "Cook X Victimization", "Married X Victimization",
                                        "Janitor X Victimization", "Paid Taxes X Victimization",
                                        "Female", "Female X Victimization", 
                                        "Youth", "Youth X Victimization", 
                                        "Tribal Member", "Tribal Member X Victimization")) %>% 
           fct_rev()) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_text(aes(label = round(estimate, 2)), vjust = -.8, size = 3) +
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2, alpha = .2) +
  labs(x = "", y = "Change: Harshness of Punishment (1-5)", subtitle = "") + 
  theme_bw()
fig9a

fig9b = regb5_2 %>% 
  tidy %>% 
  filter(!str_detect(term, "Intercept")) %>% 
  mutate(term = term %>% str_remove("Act") %>% 
           str_remove("Woman") %>% str_remove("Youth0") %>% 
           str_remove("Tribal_Member") %>% str_remove("Victimization0"),
         term = recode(term, "Victimization:Cook for Fighters" = "Cook X Victimization",
                       "Victimization:Married a Fighter" = "Married X Victimization",
                       "Victimization:Janitor at Municipality" = "Janitor X Victimization",
                       "Victimization:Paid Taxes" = "Paid Taxes X Victimization",
                       "Victimization:Female" = "Female X Victimization",
                       "Victimization:Youth" = "Youth X Victimization",
                       "Victimization:Tribal Member" = "Tribal Member X Victimization"
         )) %>% 
  mutate(term = factor(term, levels = c("Cook for Fighters", "Married a Fighter",
                                        "Janitor at Municipality", "Paid Taxes", "Victimization",
                                        "Cook X Victimization", "Married X Victimization",
                                        "Janitor X Victimization", "Paid Taxes X Victimization",
                                        "Female", "Female X Victimization", 
                                        "Youth", "Youth X Victimization", 
                                        "Tribal Member", "Tribal Member X Victimization")) %>% 
           fct_rev()) %>% 
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_text(aes(label = round(estimate, 2)), vjust = -.8, size = 3) +
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = 2, alpha = .2) +
  labs(x = "", y = "Change: Willingness to Forgive (0-1)", subtitle = "") + 
  theme_bw()
fig9b

fig9_complete = fig9a + 
  fig9b + 
  plot_annotation(title = "Effects of Respondent Victimization, Base of IS Taxpayer") & 
  theme(plot.title = element_text(hjust = 0.5))
fig9_complete
ggsave("figs/irq_fig9.pdf", fig9_complete, height = 6, width = 10)
